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I.  Introduction 


Among  the  factors  contributing  to  the  wear  and  erosion  of  gun  tubes, 
the  leakage  of  the  high  temperature  propellant  gases  past  the  moving 
projectile  occupies  an  important  role.  Presently,  there  doesn't  exist 
a rational  basis  for  the  estimation  of  the  gas  leakage  for  any  given 
gun  configuration.  This  report  describes  the  physics  of  the  leakage 
process  and  outlines  a program  for  a quantitative  determination  of  the 
leakage  flow  for  smooth  bore  gun  tubes. 

In  the  following  sections  we  describe  the  physical  phenomena,  the 
assumptions  of  the  mathematical  model,  and  give  in  an  abbreviated  form 
the  governing  equations.  We  then  discuss  the  simplifications  which  have 
to  be  made  to  enable  the  development  of  an  efficient  algorithm  which  is 
then  outlined.  The  report  closes  with  recommendations  for  implementa- 
tion of  the  outlined  work. 

II.  Description  of  the  Physical  Process 

Assume  that  a projectile  in  a smooth  bore  gun  tube  at  time  t = 0 
starts  to  accelerate  at  the  breech  end  due  to  the  pressure  rise  of  the 
gases  generated  in  the  propellant  bed  by  combustion.  The  gases  are  com- 
pressible, viscid,  turbulent  and  contain  unburnt  propellant  fragments. 

The  ignition  phenomena  leading  to  the  gas  generation  are  complex 
and  still  need  to  be  explored  in  fine  detail.  The  sequence  of  events 
is  pictured  as  follows:  first,  the  base  pad  is  ignited  by  the  primer 
gases.  The  center  core  ignites  and  begins  to  spread  the  flame.  The 
charge  itself  is  then  ignited,  the  propellant  bag  breaks  and  the  bed 
expands  to  fill  the  chamber,  thereby  transforming  the  latter  from  a 
packed  to  a fluidized  condition.  See  Figure  1.  As  the  pressure  rises 
due  to  the  burning,  the  bore  resistance  is  overcome  and  the  projectile 
begins  to  move. 

Initially  the  two  phase  mixture  consisting  of  product  gases  and 
unburnt  particles  is  moving  subsonically  with  respect  to  the  tube  wall. 
However,  it  and  the  projectile  continuously  accelerate  reaching  a velo- 

city  of  the  order  of  10  (m/s)  at  the  gun  tube  exit.  The  temperature  of 
the  gases  near  the  muzzle  is  around  2000° (K),  while  the  wall  is  initial- 
ly at  ambient  conditions.  The  gas  motion  is  highly  turbulent  and  un- 
steady. Due  to  the  presence  of  unburnt  particles  and  solid  combustion 
products,  the  ratio  of  the  specific  heats  of  the  mixture  of  particles 
and  gases  is  lower  than  that  for  the  combustion  product  gases  alone. 

Through  the  clearance  between  the  projectile  and  the  gun  tube  wall, 
which  is  of  the  order  of  0.2  mm  in  a new  tube,  the  high  pressure  gases 
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from  the  region  behind  the  projectile  leak  past  and  outrace  the  moving 

projectile.  Typically,  as  shown  by  Loder^,  for  the  105  mm  gun  tube,  at 
0.5  calibers  from  the  beginning  of  the  rifling,  which  is  also  near  the 
position  of  maximum  wear,  the  groove  due  to  the  wear  can  be  up  to  2.28 
mm  (90  mils)  measured  land  to  land.  The  severity  of  the  wear  decreases 
as  one  goes  from  the  breech  toward  the  muzzle  end  of  the  tube.  Interest- 
ingly enough,  wear  just  as  severe  has  been  reported  for  smooth  bore  guns. 

2 

For  a summary  of  our  current  understanding  of  wear  see  Ward  . 

The  nature  of  the  flow  around  the  projectile  in  the  gun  tube  will 
be  directly  affected  by  the  previous  history  of  bore  wear,  i.e.,  the 
change  in  the  area  available  for  the  passage  of  the  two  phase  flow.  The 
rate  of  leakage  will  depend  on  the  pressures  behind  and  ahead  of  the 
projectile  and  of  course  the  clearance  between  the  projectile  and  the 
gun  tube.  The  pressure  of  course  depends  on  the  time  and  the  position, 
while  the  clearance  is  position  dependent.  Hence,  the  leakage  is  time 
and  position  dependent.  That  leakage  can  be  appreciable  was  already 

3 

observed  by  Oswatitsch  , who  found  that  the  pressure  ahead  of  the  pro- 
jectile in  a gun  tube  can  be  as  high  as  10%  of  the  pressure  prevailing 
at  its  base. 

The  physical  dimensions  and  time  scales  of  the  phenomena  under  con- 
sideration are  as  follows.  The  projectile  is  assumed  to  have  a length 
to  diameter  ratio  exceeding  20.  It  accelerates  from  rest  and  at  the 
muzzle  travels  at  a low  supersonic  velocity.  Typical  transit  times 
for  a medium  to  large  caliber  weapon  are  therefore  of  the  order  of  10  - 
20  ms.  The  gap  between  the  projectile  and  the  barrel  is  of  the  order 
of  0.2  - 5%  of  the  barrel  diameter. 

While  the  pressure  behind  the  projectile  is  in  the  megapascal 
regime,  locally  encountered  temperatures  hover  around  2000  [K] . Typi- 
cal Reynolds  numbers  based  on  the  barrel  diameter  are  of  the  order  of 

108. 

The  model  adopted  here  assumes  high  pressure  gases  leaking  past  a 
moving  projectile  which  in  turn  moves  with  respect  to  a smooth  bore  gun 
tube.  The  gases  will  be  taken  to  be  premixed  and  of  a single  phase. 

In  reality  of  course,  one  has  to  deal  with  a two  phase  flow.  To  include 
some  of  the  two  phase  effects  into  a single  phase  description,  a ratio 
of  the  specific  heat  of  the  gas  will  be  used  which  is  lower  than  that 
of  a pure  gas. 


Loder,  R.K.;  Ballistic  Research  Laboratory,  Private  Communication, 

August  1977. 

^Ward,  J.R.,  "A  New  Initiative  in  Gun  Barrel  Wear  and  Erosion",  Proceed- 
ings of  the  Tri-Service  Gun  Tube  Wear  and  Erosion  Symposium,  US  Army 
Armament  Research  and  Development  Command,  Dover,  NJ,  1977. 

30swatitsch,  K.,  "Intermediate  Ballistics",  DVL  Rpt.  No.  358,  Aachen,  1964. 
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The  mixture  is  allowed  to  take  part  in  a chemical  transformation. 

A simple  decomposition  scheme  A-*B  is  adopted,  where  all  the  reactants 
are  lumped  into  A and  the  reaction  products  into  B.  The  latter  are 
assumed  to  have  attained  chemical  equilibrium.  Finite  rate  chemistry 
is  assumed  to  control  the  disappearance  of  the  reactants  with  a rate 
parameter  varying  between  1 and  0,  with  1 indicating  the  presence  of 
reactants  only  and  0 the  completion  of  the  chemical  change.  The  exo- 
thermicity  of  the  reaction  will  be  accounted  for  through  the  enthalpy 
change  of  the  gases.  Arrhenius  kinetics,  it  is  felt,  describes  the 
chemical  evolution  of  the  system  adequately  with  the  reaction  rate  con- 
stant and  the  activation  energy  taken  as  a composite  of  the  most  domi- 
nant species. 

To  adequately  model  the  flow,  one  must  adopt  a viable  turbulence 
model.  Assuming  the  validity  of  Morkovin's  hypothesis,  which  states 
that  the  effect  of  compressibility  upon  turbulence  is  negligible,  ex- 
cepting for  an  influence  on  the  mean  density  of  the  gas,  the  k-e  model 

4 

of  turbulence  of  Launder  and  Jones  , proven  for  incompressible  flow 
calculations,  will  be  used. 

The  possibility  of  relaminarization  of  the  flow,  as  it  emerges  from 
the  narrow  channel  between  the  projectile  and  the  gun  tube  wall,  exists 

and  will  be  admitted.  The  acceleration  parameter  k = (v/u2)  (du/dx) > 10  ' 

a 

for  typical  flows  in  gun  tubes  indicates  the  likelihood  of  relaminari- 
zation. Here  v is  the  kinematic  viscosity  and  u the  flow  velocity. 

This  turn  of  events  is  thought  to  come  about  because  the  turbulent 
boundary  layer  is  unable  to  follow  the  acceleration  and  with  the  fall 
of  the  Reynolds  number,  viscous  effects  assume  enough  importance  within 
the  boundary  layer  to  lead  to  the  growth  and  the  ultimate  dominance  of 
the  sublayer. 

III.  The  Proposed  Mathematical  Model 
A.  General  Observations 

Considerable  insight  and  mathematical  simplification  results 
from  an  examination  of  the  time  constants  of  the  dominant  processes  of 
the  problem.  In  order  to  establish  some  idea  of  the  time  scale  for  the 
diffusion  of  heat,  consider  the  highly  simplified  case  of  pure  diffusion 
(i.e.,  no  convection)  of  temperature  T in  one  space  dimension.  The  ap- 
propriate equation  can  be  written 


4 Launder,  B.  E.  and  Spalding,  D.  B.,  Lectures  in  Mathematical  Models 
of  Turbulence,  Academic  Press,  London,  1972. 


11 


with  boundary  conditions 


y > 0 

t .<  0 

T(y,t)  = T 

II 

o 

t 0 

T(o,t)  = 0 

and  this  equation  has  the  solution 


T = 


T 


o 


dr. 


where 


n = y/(4vt/Pr)1/2 


Here  v is  the  usual  laminar  kinematic  viscosity  and  Pr  is  the  conven- 
tionally defined  Prandtl  number.  The  solution  represents  the  transient 
growth  of  the  cold  wall  temperature  defect  into  a hot  external  environ- 
ment. If  the  scale  of  length  y in  the  definition  of  v is  taken  to  be 
the  tube  radius  and  time  scale  the  time  for  the  gas  to  pass  the  projec- 
tile. then  typically  n would  achieve  values  of  fifty  and  higher.  It 
can  readily  be  deduced  from  the  previously  given  temperature-time  solu- 
tion that  using  these  typical  values  of  n would  result  in  thermal  bound- 
ary layers  which  are  extremely  thin  c.nd  of  course  transiently  growing. 

In  the  problem  of  interest  here,  of  course,  such  thin  transiently  vary- 
ing thermal  boundary  layers  control  the  wall  heat  transfer  and  cannot  be 
neglected. 

The  fluid  flow  equations  of  motion  given  earlier  contain  time 
derivatives  of  the  density,  velocitv  and  enthalpy  of  the  gas.  If  these 
time  derivatives  were  all  negligible  compared  to  the  spatial  derivatives, 
a quasi-steady  flow  approximation  would  be  a reasonable  assumption  which 
could  be  used  to  reduce  the  computational  labor.  Order  of  magnitude 
considerations  lead  to  the  conclusion,  for  instance,  that  the  temporal 
density  gradient  is  important  and  cannot  be  neglected.  Hence, as  might 
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be  supposed  on  intuitive  grounds,  the  transient  effects  on  the  aerody- 
namic process  within  the  barrel  are  not  negligible  either.  Thus  during 
most  of  the  transit  time  the  flow  should  exhibit  significant  transient 
effects  on  the  aerodynamics  and  heat  transfer  and  this  is  substantiated 
by  (among  others)  the  calculations  of  Anderson,  et  al  (Reference  5). 

It  is  probably  self-evident  that  chemical  reactions  have  a major 
role  to  play  in  the  problem.  The  obvious  one  is  of  course  the  propel- 
lant burn  after  ignition  and  the  prospect  of  gas  phase  reaction  in  the 
gap  itself  is  feasible  and  of  great  concern.  Further,  should  the  burn- 
ing be  very  rapid  and  the  gas  be  in  or  near  chemical  equilibrium  the 
proper  determination  of  the  gap  heat  transfer  rates  requires  that  the 
change  in  chemical  equilibrium  conditions  due  to  changing  flow  environ- 
ment be  taken  into  account.  Thus, in  addition  to  the  other  previously 
mentioned  factors  we  are  led  to  consider  a (simple)  multicomponent 
reacting  flow.  For  the  leakage  problem  then,  at  least  an  axially  sym- 
metric, viscous  (turbulent),  time-dependent,  reacting  compressible  flow 
must  be  considered. 

The  presence  of  the  base  corner  gives  rise  to  very  significant 
flow  gradient  along  the  barrel,  in  addition  to  the  conventional  gradi- 
ents that  occur  in  most  boundary  layer  type  of  problems  and  which  ozeu. 
normal  to  the  wall.  Diffusion,  i.e.,  transport  effects,  upon  heat  mass 
and  momentum  in  both  the  radial  and  the  axial  direction  must  be  allowed 
for  near  the  gap  inlet.  The  pressure  field  at  the  gap  entrance  and  exit 
is  also  not  readily  approximated  or  neglected,  leading  one  to  consider 
the  full  Navier-Stokes  equations  as  being  required  to  describe  the  flow. 

B.  Governing  Equations 

The  flow  under  consideration  is  a turbulent  chemically  react- 
ing multi-component  mixture  with  heat  and  mass  transport.  The  govern- 
ing system  of  partial  differential  equations  describing  this  process  is 
based  on  the  conversion  laws  of  mass,  momentum,  energy  and  chemical 

species6.  For  simplicity  these  equations  are  expressed  in  vector  nota- 
tion below  and  all  quantities  are  nondimensional . Velocities  are  nor- 
malized by  Up,  density  by  Pp,  enthalpy  by  hp,  temperature  by  Tp,  molecu- 
lar weights  by  Wp,  pressure  by  Pp  = poRgTpZp,(Zp  = 103/'Wp) , dynamic 


^Anderson,  L.  W. , Bartlett,  E.  P.,  Dahm,  T.  J.  and  Kendall,  R.  M. , 
"Numerical  Solution  of  the  Non-Steady  Boundary  Layer  Equations  with 
Application  to  Convective  Heat  Transfer  in  Guns",  Aerotherm  Final  Rpt . 
No.  70-22,  October  1970. 

6Bird,  R.  B.,  Stewart,  W.  E.  and  Lightfoot,  E.  N.,  Transport  Phenomena, 
Wiley,  New  York,  1960. 
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viscosity  by  y , radiation  energy  flux  by  q , and  time  by  (L/U  ) where 

L is  the  reference  length.  Coupling  between  concentration  and  thermal 
gradients  (Soret  and  Dufour  effects),  pressure  gradient  diffusion,  body 
forces  and  bulk  viscosity  are  all  assumed  to  be  negligible.  In  addition, 
Fick's  law  is  presumed  valid  which  implies  equal  binary  diffusion  coef- 
ficients for  each  pair  of  species  in  the  mixture  (see,  e.g..  Reference  7). 
The  resulting  set  of  ensemble-averaged  equations  is 

Continuity 


3p 

at 


-V  • (pu) 


(i) 


Conservation  of  Gas  Phase  Species 


3(pmi)  _ 1 

:pjr  = -V*  (punu ) + — V-d^Vnu)  ♦ r.  + Es.  ^ 


J 


i.l 


(2) 


Conservation  of  Momentum 


, .,.(pSa,  . _s_  w . i_  ,.(2^5)  - |i_  n„eff(,.a)] 


at 


»DUD 


(3) 


Williams,  F.  A.,  Combustion  Theory.  Addison-Weslev.  Reading.  MA . 
1965. 


I 


Conservation  of  Energy 


*1PHI=  _y.(paH)  ♦ jL 


v*crhvn) 


* Re  h“  ’-[("eff  - ’h)’  (U2U)j 

+ ^ V-qR-^r.hj 


(4) 


p is  the  density,  u the  velocity  vector,  m.  is  the  i^*1  component  gas,  f. 
component  of  other  phase  (liquid  or  solid),  p is  the  pressure,  u J 


stagnation  enthalpy.  Re  the  Reynolds  number  based  on 


H the 
reference  length 


and  reference  velocity,  F^,  Tm,  r^,  are  the  turbulent  exchange  coeffici- 


ents of  heat,  gas  and  the  other  phase,  e is  the  mean  flow  rate  of 
strain  tensor.  A turbulence  model  is  employed  to  define  the  effective 
velocity,  y0££-  Similarly,  a chemistry  model  is  employed  to  specify  the 


production  rate  r^  for  the  chemistry,  a phase  change  model  could  be  em- 
ployed for  the  source  terms  R.  and  s.  .,  and  a radiation  transport  model 
again  could  be  employed  to  specify  These  terms  are  included  for 

completeness  but  will  not  be  used  at  least  during  the  initial  phase  of 
the  investigation,  apart  of  course  from  the  turbulence  terms. 


To  account  for  the  turbulent  behavior  in  the  solution  of  the  ensemble- 
averaged  Navier-Stokes  equations,  a turbulence  model  is  introduced  to 
define  an  effective  viscosity.  A review  of  turbulence  models  is^avail- 
able  in  the  literature  (see,  e.g.,  References  4 and  8).  Prandtly  was 
perhaps  the  first  to  introduce  a turbulence  model  when  he  postulated 
that  the  time-averaged  shear  stress  and  the  time-averaged  velocity  gradi- 
ent are  proportional  as  in  laminar  flow,  and  that  the  length  scale  (the 
so-called  mixing  length)  which  enters  the  relationship  is  proportional 
to  the  turbulent  shear  region  thickness.  A disadvantage  of  the  mixing 
length  model  is  that  it  is  an  equilibrium  model  (i.e.,  turbulence  is 
assumed  to  be  produced  and  dissipated  locally)  and  it  requires  an  ad  hoc 
mixing  length  distribution.  Some  of  the  shortcomings  of  the  mixing  length 
model  have  been  overcome  for  many  cases  of  interest  by  the  introduction 
of  various  multiequation  transport  models  of  turbulence. 


8Harlow,  F.H.,  ed. , "Turbulence  Transport  Modeling",  AIAA  Selected  Reprint 
Series,  Vol.  XIV,  1973. 

9Prandtl,  L.,  "Bericht  Uber  Untersuchungen  Zur  Ausgebildeten  Turbulenz," 
ZAMM,  Vol.  5,  1925,  p.  136. 
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Many  of  the  two-equation  turbulence  models  employ  the  Prandtl- 
Kolmogorov  formula  for  specification  of  the  turbulent  viscosity,  u^, 

57  = C pk  l (6) 

Re  y t 

where  k is  the  turbulence  kinetic  energy  and  £ is  a length  scale  of 
the  turbulence.  This  relation  follows  from  dimensional  arguments  for 
turbulent  flow  described  by  the  two  parameters,  k and  i . A major 
advantage  of  a two-equation  turbulence  model  compared  to  a mixing  length 
model  is  the  fact  that  the  length  scale  (as  well  as  the  turbulence  ki- 
netic energy)  is  determined  from  transport  equations,  whereas  in  mixing 
length  models  the  length  scale  is  determined  from  an  ad  hoc  algebraic 
expression.  Successful  use  of  the  mixing  length  model  requires  an  a 
priori  specification  of  the  turbulence  length  scale.  While  a realistic 
assumption  can  be  made  for  certain  flows  such  as  shear  layers,  in  many 
flows  of  interest  (such  as  internal  recirculating  flows)  the  choice  of 
a proper  turbulence  length  scale  is  not  obvious.  Since  the  turbulence 
length  scale  emerges  from  the  solution  in  a two-equation  model,  these 
models  are  more  likely  to  give  accurate  predictions  over  a wide  range 
of  geometric  and  flow  conditions  with  the  same  empirical  constants. 

It  should  be  noted  that  the  two-equation  models  employ  the  eddy- 
viscosity  formulation  for  the  Reynolds  stresses  as  in  the  mixing  length 
model , i . e. , , a,. 


ou.  'u'.  = — 

1 j Re 


Hence,  this  formulation  still  suffers  from  the  physical  shortcoming  that 
there  is  zero  Reynolds  stress  wherever  the  velocity  gradient  is  zero. 

In  addition,  the  eddy  viscosity  formulation  is  isotropic  which  may  be 
incorrect  in  many  three-dimensional  and  swirling  flows.  However,  for 
practical  calculations  of  complex  turbulent  internal  flows  there  are  no 
other  available  transport  models  which  are  as  suitable  or  even  as  rela- 
tively well  developed. 

Various  forms  of  the  two-equation  model  of  turbulence  have  been  pro- 
posed since  Kolmogorov1  first  introduced  the  concept  in  1942.  Most 
investigators  have  chosen  the  kinetic  energy  of  turbulence,  k,  as  their 
first  variable.  However,  there  is  a wide  diversity  of  choice  as  to  the 
second  variable  to  be  used.  In  general,  each  investigator  chose  a second 
variable  which  he  felt  was  appropriate  to  the  physical  description  of 
turbulence.  For  instance,  Kolmogorov  chose  as  his  second  variable  a 


Kolmogorov,  A.N.,  "Equations  of  Turbulent  Motion  of  an  Incompressible 
Turbulent  Fluid",  izv>  Akad.  Nauk,  SSR  Ser.  Phys.  VI,  No.  1-2,  56  1942. 
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quantity  which  is  proportional  to  the  mean  frequency  of  the  most  ener- 

11  12 

getic  motions,  while  Spalding  and  Saffman  chose  a quantity  that  rep- 
resented the  time-averaged  square  of  the  vorticity  fluctuations.  Another 
commonly  chosen  second  variable  has  been  the  turbulence  kinetic  energy 
dissipation  rate,  e,  which  was  selected  for  this  investigation.  An  ad- 
vantage of  using  the  dissipation  rate  of  turbulence  kinetic  energy  for 
the  second  equation  is  that  the  dissipation  rate  appears  directly  in 
the  turbulence  kinetic  energy  equation,  and  that  an  equation  for  e can 
be  readily  developed.  Since  derivations  of  the  equations  for  k and  e 
are  lengthy  and  have  previously  been  presented  in  the  open  literature, 
these  derivations  are  not  repeated  in  the  present  report.  The  appropri- 
ate transport  equations  for  turbulent  kinetic  energy  and  energy  dissipa- 
tion rate,  valid  at  high  Reynolds  numbers,  are  taken  directly  from  Launder 
13 

and  Spalding  . 

The  turbulence  kinetic  energy  equation  in  vector  notation  is 


= -V*  (pfik)  + 


(8) 


where  the  first  two  terms  on  the  right-hand  side  represent  convection 
and  diffusion  of  turbulence  kinetic  energy,  respectively;  the  third  and 
fourth  terms  represent  the  generation  (due  to  shear  forces)  and  dissipa- 
tion of  turbulence  kinetic  energy  respectively.  The  equation  for  the 
dissipation  rate  of  turbulence  kinetic  energy  is 


3(pe)  = 
3t 


= -V*(pue)  + 
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where  the  function  of  the  terms  is  analogous  to  thatin  the  turbulence 
kinetic  energy  equation.  In  Equations  (8)  and  (9),  e is  the  mean  flow 
rate  of  strain  tensor.  Using  dimensional  arguments  the  Prandtl- Kolmo- 
gorov formula.  Equation  (6),  may  be  written  as 


(10) 


which  implies  that  a turbulence  length  scale  or  "mixing  length"  may  be 
defined  as 


^Spalding,  D.  B.,  "The  Prediction  of  Two-Dimensional,  Steady  Turbulent 
Flows",  Imperial  College,  Heat  Transfer  Section  Report  EF/TN/A/16,  1969 

12 

Saffman,  P.  G.,  "A  Model  for  Inhomogeneous  Turbulent  Flow",  Proc.  Roy. 
Soc . , Vol . A317,  1970,  p.  417. 

13 

Launder,  B.  E.  and  Spalding,  D.  B.,  "The  Numerical  Computation  of  Tur- 
bulent Flows",  Computer  Methods  in  Applied  Mechanics  and  Engineering, 
Vol.  3,  1974,  p.  269. 
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The  constants  appearing  in  Equations  (8)  through  (10)  have  been 

14 

evaluated  by  Jones  and  Launder  and  the  values  proposed  are 


IV. 


Cp  = 0.09 
C = 1.55 
C2  = 2.0 

°k  - 1-° 
a =1.3 


e 

Discussion  of  the  Numerical  Treatment 
A.  General  Observations 


(12) 


Having  concluded  in  the  previous  section  that  a rigorous  treat- 
ment of  the  physics  of  the  leakage  problem  demands  that  the  full  Navier- 
Stokes  equations  be  used,  it  now  becomes  necessary  to  discuss  the  prob- 
lems involved  in  obtaining  numerical  solutions  to  these  equations.  The 
full  Navier-Stokes  equations  do  have  the  attractive  features  that  at 
least  in  a formal  sense  their  use  minimizes  the  physical  approximations 
required.  However,  both  physical  and  numerical  treatments  are  required 
to  model  on  the  one  hand  the  turbulent  transport  and  on  the  other  to 
deal  with  the  coupled  nonlinear  nature  of  the  Navier-Stokes  equations. 
Here  the  numerical  problems  which  arise  when  the  Navier-Stokes  equations 
are  solved  under  the  conditions  of  the  leakage  problem  are  discussed. 

The  chosen  numerical  treatment  is  a critical  factor  in  determin- 
ing the  economics  or  perhaps  even  the  feasibility  of  performing  the  de- 
sired calculations.  Use  of  a numerical  algorithm  not  well  suited  to  the 
problem  at  hand  could  result  in  a prohibitive  run  time  to  generate  solu- 
tions to  the  required  degree  of  accuracy.  Consideration  of  the  need  to 
accurately  define  the  very  thin  boundary  layers,  both  thermal  and  momen- 
tum, leads  to  proposing  the  use  of  an  implicit  finite  difference  algo- 
rithm. Explicit  and/or  stability  restricted  schemes  are  usually  forced 
to  take  extremely  small  time  steps  as  a result  of  the  refined  spatial 
mesh.  This  unduly  small  time  step  is  much  smaller  than  would  be  required 
for  transient  accuracy.  The  implicit  schemes  usually  do  not  suffer  from 
this  particular  problem  and  can  be  very  much  more  efficient  for  boundary 
region  type  of  problems.  Implicit  schemes,  however,  have  not  been  widely 
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used  for  transient  mixed  convection-diffusion  problems  and  the  possible 
problems  which  might  arise  in  this  area  are  discussed  in  detail.  Fur- 
ther, the  choice  of  an  implicit  scheme  leads  to  the  necessity  of  consid- 
ering in  detail  equation  linearization  and  coupling,  and  this  will  be 
done.  A Taylor  series  expansion  in  time  leads  to  a very  powerful  lin- 
earization and  this  concept  is  recommended.  With  an  explicit  scheme 
the  choice  of  dependent  variables  can  virtually  eliminate  the  problem 
of  linearization  and  coupling  which  is  why  in  the  literature  this  prob- 
lem is  often  not,  as  yet,  well  discussed. 

Common  to  any  algorithm  is  the  coordinate  system  and  the  boundary 
conditions.  The  choice  of  coordinate  system  and  the  subsequently  used 
mesh  can  be  a critical  factor  in  the  overall  algorithm  accuracy.  Conse- 
quently both  the  coordinate  system  and  the  mesh  will  be  discussed  in 
detail.  On  the  basis  of  this  discussion  it  will  be  suggested  that  a 
body  fitted  coordinate  system  based  upon  the  projectile  would  be  both 
convenient  and  accurate. 

B.  Transient  Accuracy,  Stability  Limits 

As  a general  observation  great  care  is  required  to  obtain  accept- 
able transient  accuracy  for  long  time  integration  of  fluid  flow  equations 
in  which  the  convective  terms  are  significant.  The  problem  is  discussed 

for  instance  by  Morton^5  and  Roache^.  In  the  more  widely  used  finite 
difference  schemes  two  types  of  transient  errors  are  commonly  encountered, 
damping  and  dispersion.  Considering  the  pure  convection  of  a wave,  damp- 
ing would  be  reflected  in  the  change  with  time  of  the  wave  amplitude 
while  dispersion  would  be  reflected  in  the  relative  phase  change  with 
time,  i.e.,  the  change  with  time  of  the  wave  propagation  velocity.  Time 
centering  of  the  spatial  derivatives  which  would  require  taking  all  spa- 
tial derivatives  at  the  midpoint  of  the  time  step,  can  be  employed  to 
remove  the  numerical  damping.  However,  dispersion  can  be  much  more  dif- 
ficult to  treat  adequately  in  order  to  achieve  a desired  level  of  accu- 
racy due  to  the  presence  of  phase  errors.  To  date  the  more  widely  used 
schemes  with  good  properties  with  regards  to  dispersion  have  been  stabil- 
ity restricted  in  the  manner  described  by  Courant,  Friedrichs  and  Lewy1^. 
In  this  sense  the  maximum  time  step  permissible  in  the  simple  compres- 
sible one-dimensional  fluid  flow  problem  is  found  to  be  restricted  such 
that 


^Morton,  K.  W. , "Stability  and  Convergence  in  FI  id  Flow  Problems", 

Proc.  Royal  Soc.,  Vol.  A323,  No.  1553,  June  1971,  pp.  237-253. 

^Roache,  P.  J.,  "Computational  Fluid  Dynamics",  Hermosa  Publishers,  1972. 

^Courant,  R.,  Friedrichs,  K.  0.  and  Lewy,  H.,  "Uber  die  Partiellen 
Differenzengleichungen  der  Mathematischen  Physik",  Mathematische 
Annalen,  Vol.  100,  1928,  pp.  32-74. 


At  < Ax/(|u|  + a)  where  u is  the  convection  velocity  and  a is  the 

local  speed  of  sound  and  Ax  is  the  spatial  mesh  increment.  Experience 
indicates  that  this  type  of  restriction  appears  even  in  the  multidimen- 
sional analogues.  Since  locally  refined  spatial  meshes  are  required  in 
the  present  problem  to  define  the  boundary  layers  and  also  as  a result 
of  the  low  local  Mach  numbers  expected  (high  temperatures),  the  Courant- 
Friedrich-Lewy  (CFL)  requirement  could  be  painfully  restrictive.  In  such 
a situation  the  implication  is  that  the  physical  processes  are  changing 
with  time  more  slowly  than  the  CFL  limit  and  that  the  time  truncation 
errors  would  be  acceptable  if  a time  step  greater  than  the  CFL  limit  were 
used  and  the  numerical  scheme  could  accept  such  a large  time  step.  The 
counter  argument  sometimes  given  is  that  indeed  the  CFL  limit  represents 
the  time  scale  of  propagation  of  physical  information;  but  this  is  clearly 
not  generally  correct  since  incompressible  flow  simulations  of  time- 
dependent  low  Mach  number  compressible  flows  are  normally  very  good  and 
in  these  analyses  the  speed  of  sound  is  infinite  with  a resulting  allow- 
able CFL  limited  time  step  of  zero.  Having  said  this,  it  should  be 
pointed  out  that  in  many  physical  processes  the  relevant  time  scale  of 
information  propagation  is  indeed  the  local  speed  of  sound  and  as  such 
the  CFL  limit  would  not  be  restrictive  (unless  only  the  steady  solution 
were  being  sought) . For  the  particular  problem  under  consideration  un- 
fortunately the  CFL  limit  would  be  unduly  restrictive  certainly  during 
the  initial  phases  of  the  motion  and  a scheme  not  restricted  by  the  CFL 
limit  would  be  preferable.  This  leads  to  the  previously  mentioned  recom- 
mendation of  an  implicit  difference  scheme. 

Turning  to  these  unconditionally  stable  schemes  Morton15  shows  the 
relatively  poor  dispersion  characteristics  of  a Crank -Nicolson  implicit 
scheme  applied  to  the  pure  convection  problem.  Although  not  shown  by 
Morton,  the  dispersion  at  a given  wavelength  increases  progressively  as 
the  CFL  limit  is  exceeded.  Thus  it  appears  that  for  pure  convection 
with  a Crank-Nicolson  implicit  scheme  taking  large  time  steps,  a wave 
described  by  less  than  ten  or  so  grid  points  would  suffer  unacceptable 
dispersion  after  some  period  of  time.  Thus  accurate  long  time  calcula- 
tions of  short  wave  length  components  of  the  motion  would  require  a very 
refined  spatial  mesh.  It  would  then  seem  that  something  of  a dilemma 
exists.  Schemes  with  good  dispersion  properties  require  many  time  steps 
to  describe  the  motion  due  to  the  CFL  restriction  (the  start  of  the 
motion  seems  particularly  difficult  with  a CFL  restricted  scheme) . 
Schemes  not  restricted  by  the  CFL  limit  do  not  require  as  many  time 
steps  but  do  require  good  spatial  resolution  in  order  to  have  acceptable 
short  wave  length  dispersion.  Some  observations  are  pertinent  at  this 
point  however.  First  is  that  good  spatial  resolution  is  required  in  any 
event  to  define  the  boundary  layers  and  shear  layers  in  the  particular 
problem  under  consideration.  Next  in  the  particular  problem  that  is 
being  considered,  the  time  dependency  is  continuously  imparted  through 
the  initial  and  boundary  conditions.  As  a result  the  problem  is  more 
one  of  predicting  the  diffusion  of  a forced  transient  than  trying  to 
determine  the  phase  lag  of  a wave  in  a slightly  sheared  flow  after  a 


long  propagation  time.  In  addition,  as  Morton  points  out  the  better 
explicit  schemes  with  good  dispersion  properties  have  certain  other  un- 
desirable properties  which  require  special  treatment  and  which  in  a com- 
plicated multidimensional  problem  might  prove  difficult  to  implement. 
Lastly,  apparently  there  have  not  been  any  significant  attempts  made  to 
improve  the  dispersion  characteristics  of  the  implicit  schemes,  and 
should  the  potential  gains  warrant  the  labor,  it  would  seem  very  worth- 
while to  investigate  a modification  to  the  implicit  scheme  similar  to 
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the  one  developed  by  Fromm  for  the  Lax-Wendroff  scheme  or  apply  the 

19 

ideas  of  Boris  and  Book  which,  although  developed  for  treating  large 
gradients  in  the  flow,  do  reduce  dispersive  errors. 

With  the  foregoing  in  mind  it  would  be  reasonable  to  proceed  with 
an  implicit  scheme  in  view  of  the  restrictions  that  might  otherwise 
result  in  the  use  of  the  stability  restricted  explicit  schemes  with  both 
low  Mach  numbers  and  refined  spatial  meshes.  Dispersion  might  be  regard- 
ed as  a potential  problem  for  the  implicit  scheme  and  therefore  both  mesh 
refinement  studies  and  model  problem  evaluation  must  be  carried  out,  both 
for  the  actual  cases  involved  and  for  a model  system  more  representative 
of  the  mixed  convection-diffusion  problem  under  consideration.  Should 
critical  flow  features  be  found  subject  to  undue  dispersion,  considera- 
tion should  be  given  to  either  increased  spatial  resolution  or  to  devel- 
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oping  a correction  for  the  implicit  scheme  following  Fromm  or  Boris 
and  Book 

As  was  mentioned  earlier,  time  centering  or  "Crank-Nicolson"  averag- 
ing of  the  spatial  derivatives  produces  a scheme  which  has  no  damping  in 
the  simple  problem  and  is  of  formal  second  order  accuracy  in  time.  For 
the  more  complicated  systems  of  equations  which  govern  the  fluid  flow 
great  care  is  required,  for  example,  with  the  treatment  of  nonlinearities 
and  cross  derivatives  to  obtain  formal  second  order  temporal  accuracy. 
When  this  is  compounded  with  the  difficulties  of  using  a highly  non- 
linear turbulence  and  chemistry  model  coupled  implicitly  with  the  depen- 
dent variables  in  a consistent  second  order  time-dependent  manner,  a 
great  deal  of  work  is  clearly  required  to  achieve  formal  second  order 
transient  accuracy  in  these  non-model  problems.  At  this  early  stage  of 
the  overall  problem  development  it  is  doubtful  if  a great  deal  of  extra 
effort  required  to  achieve  formal  second  order  transient  accuracy  is 
warranted.  Much,  of  course,  can  and  should  be  done  to  improve  the  tran- 
sient accuracy.  After  the  fact  and  as  part  of  the  investigation  of  the 
dispersion  problem  an  evaluation  of  the  benefit  of  more  complete  second 
order  accurate  transient  representation  can  be  made. 
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C.  Equation  Coupling  and  Linearization  Problems 

In  the  previous  section  it  was  argued  that  a conditionally  sta- 
ble explicit  scheme  would  require  a very  large  number  of  time  steps  to 
integrate  over  the  projectile  transit  time  in  view  of  the  initial  low 
Mach  number  and  locally  refined  mesh  used  in  this  mixed  convection-dif- 
fusion problem.  This  observation  led  to  the  proposed  use  of  an  implicit 
scheme  and  once  an  implicit  scheme  is  envisaged,  equation  coupling  and 
linearization  must  be  considered.  As  mentioned  earlier,  with  an  expli- 
cit scheme  the  choice  of  dependent  variables  can  mitigate  greatly  this 
particular  problem  such  that  it  really  has  come  to  be  thought  of  as  pri- 
marily a feature  of  implicit  methods. 

The  linearization  problem  leads  naturally  into  the  coupling 
problem  so  it  is  discussed  here  first  of  all.  Both  problems  are  reviewed 
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in  detail  by  McDonald  and  Briley  and  Briley  and  McDonald  . It  is 
argued  by  these  authors  that  in  order  to  get  the  most  accuracy  out  of 
a given  grid,  the  errors  arising  from  representing  nonlinear  terms  by 
linear  combinations  of  terms  at  the  unknown  time  level  should  be  less 
than  or  equal  to  the  discretization  errors.  If  the  linearization  errors 
are  greater  than  the  temporal  or  spatial  discretization  errors  then  clear- 
ly the  accuracy  of  the  differencing  is  being  squandered  and  iteration  or 
some  form  of  linearization  improvement  is  called  for.  As  a first  option, 
iteration  to  reduce  the  linearization  error  is  not'  to  be  recommended, 
since  iteration  only  improves  the  linearization;  yet  computationally 
usually  it  costs  as  much  for  one  iteration  as  it  does  to  march  one  time 
step.  Reducing  the  time  step  would  be  preferable  to  iteration  since 
cutting  back  on  the  time  step  would  improve  both  the  transient  accuracy 
and  reduce  the  linearization  error.  To  obtain  a linearization  which  in- 
troduces errors  of  at  most  the  same  order  as  the  temporal  differencing, 
a Taylor  series  expansion  about  the  known  time  level  can  be  performed, 
and  this  is  the  approach  adopted  in  References  22  and  23.  The  process 
results  in,  for  instance,  expanding  a uv  term  about  the  known  time  level 
n to  obtain  a second  order  accurate  linear  expression  at  the  unknown  time 
level  n+1. 
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and  immediately  a coupled  system  is  obtained,  since  neither  u nor  v can 
be  left  at  the  old  time  level  in  order  to  generate  an  uncoupled  system 
(i.e.,  one  scalar  equation  for  each  dependent  variable  containing  only 
that  particular  dependent  variable  at  the  unknown  time  level)  without 
introducing  a first  order  time  truncation  error.  Thus  the  formal  lin- 
earization process  and  consideration  of  the  resulting  errors  indicates 
that  from  an  accuracy  point  of  view  a coupled  system  of  equations  at 
the  new  time  level  must  be  treated.  This  can  be  efficiently  accom- 
plished by  use  of  a block  implicit  scheme,  particularly  if  the  blocks 
are  relatively  small  (e.g.,  equal  to  the  number  of  governing  partial 
differential  equations,  at  most). 


As  a related  topic  the  sequential  solution  of  ad  hoc  uncoupled  sys- 
tems is  sometimes  advocated.  This  approach  has  merit  in  cases  where  a 
weak  coupling  between  the  governing  equations  exists  and  can  be  exploited. 
Blottner“  has  shown  that  in  the  case  of  the  steady  boundary  layer  equa- 
tions such  an  ad  hoc  uncoupling  can  seriously  degrade  the  solution  accu- 
racy and  the  sometimes  more  than  ten  circuits  round  the  uncoupled  system 
is  required  to  achieve  the  accuracy  level  obtained  in  one  pass  with  the 
coupled  scheme.  Consequently  in  spite  of  the  increased  labor  involved 
in  solving  the  block  system  a substantial  net  saving  is  obtained  vis  a 
vis  the  sequential  iterative  approach.  In  the  present  problem  the  coup- 
ling effects  between  equations  can  be  much  stronger  than  are  observed 
in  the  conventional  boundary  layer  equations  and  so  the  advantage  of 
the  block  system  approach  would  probably  be  considerable. 

D.  The  Coordinate  System  and  Related  Topics 


The  projectile  in  the  barrel  has  three  distinct  geometric 
regions  of  interest.  Referring  to  Figure  2 a schematic  is  given  which 
illustrates  the  three  regions.  Region  I is  the  base  region  of  the  pro- 
jectile while  region  II  is  the  gap  itself,  shown  much  enlarged  in  the 
schematic.  Region  III  is  the  projectile  forebody  and  this  region  is 
important  in  the  present  problem  mainly  to  set  the  exit  conditions  for 
the  leakage  flow  from  region  II.  Two  convenient  coordinate  systems  for 
regions  I and  II  come  immediately  to  mind.  The  first  of  these  is  the 
use  of  an  axisymmetric  variant  of  the  blunt  body  conformal  coordinate 
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of  Davis  (see  Ghia  and  Davis  ).  In  this  two-dimensional  case  Ghia  and 
Davis  constructed  a coordinate  system  for  the  region  near  the  front  of 
a blunt  body  in  a uniform  stream  by  a conformal  transformation  scheme. 
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By  use  of  an  image  body  the  barrel  itself  could  be  introduced  and 
regions  I and  II  described  in  a conformal  coordinate  system.  At  some 
point  down  the  gap  this  coordinate  system  would  have  to  be  joined  to  a 
suitable  forebody  coordinate  system,  but  this  would  not  seem  to  be  an 
insurmountable  problem. 

The  alternative  and  recommended  scheme  for  region  I and  II  is  to 
adopt  a simple  Cartesian  mesh  stretched  in  both  r and  z by  an  analytic 
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stretching  function  such  as  that  described  by  Roberts^  . The  analytic 
stretching  function  could  then  be  used  to  place  a relatively  dense  mesh 
in  the  vicinity  of  the  corner.  As  with  the  conformal  transformation 
scheme  the  forebody  coordinates  would  have  to  be  adjoined  to  the  gap 
coordinate.  In  view  of  the  prior  success  which  investigators  have  had 
with  both  aforementioned  coordinates  there  is  no  telling  argument  for 
or  against  either  system.  The  stretched  Cartesian  mesh  is  simpler  to 
implement  and  as  a consequence  is  the  recommended  approach  at  least  for 
the  initial  stages  of  this  work. 

Joining  of  either  coordinate  systems  to  the  forebody  region  III 
could  be  performed  discretely  with  one-sided  differencing  used  on  either 
side  of  the  join.  Thus  discontinuous  derivatives  of  the  coordinates 
could  be  permitted  at  the  join  line.  While  feasible,  and  indeed  such 
discretely  joined  differencing  meshes  have  been  used  in  other  investi- 
gations, in  the  present  problem  this  join  problem  can  be  moved  to  the 
body  nose  by  normalizing  the  radial  mesh  by  the  gap  height  h(z).  The 
mesh  can  then  be  continuous  from  region  I into  region  III.  Difficulty 
now  occurs  at  the  body  nose  where  3h/3z  is  discontinuous.  The  nose 
problem  can  now  be  treated  by  one-sided  differencing  of  3h/3z  and  the 
flow  variables  across  the  discontinuity.  The  principal  attribute  of  the 
h(z)  normalization  is  that  it  provides  a very  simple  method  of  treating 
arbitrary  forebody  geometries  at  the  same  time  as  conveniently  blending 
region  II  and  region  III  meshes.  The  singularity  of  3h/3z  at  the  body 
nose  is  a problem, however,  that  may  require  special  treatment.  The 
corner  region  mesh  and  its  blend  into  region  I is  not  a problem  at 
least  insofar  as  h(z)  is  concerned  since  in  region  I h(z)  can  be  fixed 
at  the  corner  value,  h(zc). 

The  problem  is  very  simply  formulated  in  a projectile  fixed  coordi- 
nate system  with  a moving  wall  to  represent  the  barrel.  A schematic  of 
this  system  is  shown  in  Figure  3a  with  the  upstream  and  downstream  com- 
putational boundary  denoted  by  a chain  line.  This  particular  coordinate 


26 

Roberts,  G.  0.,  "Computational  Meshes  for  Boundary  Layer  Problems", 
Proceedings  of  the  Second  International  Conference  on  Numerical  Methods 
in  Fluid  Dynamics,  Springer-Verlag,  New  York,  p.  171,  1971. 


24 


4 


m 


system  has  the  attribute  that  definition  is  retained  in  the  gap  region 
without  difficulty  as  the  flow  develops  in  time.  The  alternative  and  for 
the  time  being,  not  recommended  moving  coordinate  system  could  be  intro- 
duced whereby  the  distance  from  the  breech  to  the  projectile  base,  £(t), 
could  be  used  to  normalize  the  axial  distance  as  is  shown  in  Figure  3b. 
The  new  axial  coordinate  would  be  related  to  the  physical  distance  by 
z(z,t)  = z/£  (t ) with  t ft)  describing  the  position  of  the  projectile. 

This  time  stretched  coordinate  reduces  the  axial  definition  in  the  gap 
as  the  solution  develops  in  time  and  also  continues  to  define  regions 
of  the  flow  where,  from  the  gap  leakage  problem  point  of  view,  interest 
may  no  longer  exist,  or  the  flow  may  be  adequately  described  by  simpler 
concepts  such  as  a boundary  layer  and  core  flow  analysis. 

In  order  to  maximize  the  definition  of  the  gap  flow,  a coordinate 
system  attached  only  to  the  projectile  base  and  not  tied  to  both  the 
projectile  and  the  breech  is  recommended  at  the  present  time.  Should 
the  interaction  between  the  gap  and  the  breech  conditions  prove  ulti- 
mately to  be  of  critical  importance,  it  is  a straightforward  matter  to 
revert  to  the  expanding  mesh  tied  both  to  the  projectile  and  the  breech. 

Corners  in  themselves  present  a major  problem  to  finite  difference 
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schemes.  This  problem  is  discussed  for  instance  by  Mitchell  and  by 
28  16 

Whiteman  . Roache  gives  a morphology  of  schemes  for  treating  the 
corner  singularity  with  finite  difference  procedures.  Of  particular 
note  is  the  observation  by  both  Mitchell  and  Whiteman  that  while  finite 
difference  solutions  lose  accuracy  near  re-entrant  corners,  mesh  refine- 
ment reduces  the  zone  of  contamination.  Whiteman  was  even  able  to  prove 
uniform  convergence  for  a family  of  finite  difference  solutions  for  a 
simple  slit  type  problem.  An  alternative  to  mesh  refinement  is  to  imbed 
a finite  series  solution  within  the  difference  solution  to  represent  the 
singularity.  This  method  is  also  attractive  and  has  been  used  with  suc- 
cess in  some  problems.  In  view  of  the  slow  rate  of  convergence  of  the 
refined  mesh  finite  difference  solution,  the  imbedding  technique  should 
be  investigated  further  to  improve,  if  possible,  its  generality  and  ease 
of  implementation.  For  the  time  being,  however,  the  locally  refined  fi- 
nite difference  approach  can  be  continued.  Whiteman  also  advocated  com- 
parison of  finite  element  with  finite  difference  solution  to  see  if 
there  was  a relative  advantage  for  this  type  of  singularity  and  to  date 
this  comparison  has  apparently  not  been  performed,  but  would  be  worth- 
while. 
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E.  Boundary  Conditions  Implementation 

It  has  been  plausibly  argued,  and  will  be  adopted  here,  that 
the  time  history  of  the  upstream  boundary  condition  can  be  prescribed 
adequately  for  gap  leakage  purposes  from  one-dimension  propulsion  codes 
already  in  existence  (see  for  instance  Anderson  et  al.)^  Thus  for  the 
initial  computation  of  gap  leakage  it  is  reasonable  to  assume  that  the 
time  varying  wall  velocity  will  be  given.  No  slip  conditions  will  then 
be  applied  on  the  solid  surface  boundaries.  The  surface  temperatures 
should  in  principle  be  determined  in  conjuction  with  a metal  conduction 
program.  In  the  initial  studies  the  complexity  of  performing  a metal 
conduction  calculation  at  each  time  step  can  be  avoided  and  the  surface 
temperature  and  heat  transfer  bounded  by  the  limits  of  adiabatic  wall 
(zero  thermal  gradient  at  the  wall)  and  surface  temperature  fixed  at 
time  zero  value.  Both  of  these  limiting  cases  are  readily  implemented 
and  should  be  evaluated.  The  wall  pressure  poses  something  of  a prob- 
lem and  here  it  is  suggested  that  the  conventional  steady  boundary  lay- 
er approximation  of  zero  normal  gradient  of  pressure  should  be  adequate 
in  the  present  problem,  only  in  this  instance  the  sum  of  the  pressure 
gradient  and  the  acceleration  terms  are  equated  to  zero. 
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With  the  limited  number  of  grid  points  available  it  is  desir- 
able to  place  the  downstream  boundary  as  close  to  the  projectile  as  pos- 
sible and  retain  gap  definition.  Three  strategies  then  are  open.  First, 
place  the  boundary  in  the  quiescent  air  outside  the  weapon  with  a 
highly  stretched  axial  mesh  and  apply  undisturbed  flow  as  boundary  con- 
ditions. Secondly,  compute  the  flow  ahead  of  the  projectile  by,  say,  an 
existing  inviscid  irrotational  flow  procedure  and  impose  boundary  condi- 
tions for  the  Navier-Stokes  analysis  close  to  the  projectile  nose  based 
on  these  inviscid  calculations.  This  leaves  open  the  question,  presum- 
ably answered  by  the  author  of  the  inviscid  procedure,  of  the  boundary 
conditions  for  the  inviscid  code.  Here  of  course  transformations  and 
the  use  of  a potential  function  ease  the  inviscid  flow  calculation  labor 
considerably  so  it  is  probably  feasible  in  this  context  to  use  quiescent 
air  outside  the  breech  as  the  inviscid  flow  downsteam  conditions.  Thirdly, 
apply  weak  exit  boundary  conditions  based  upon  linear  extrapolation  from 
interior  grid  points.  Here  the  terminology  "weak"  refers  to  the  con- 
straint of  the  boundary  condition  upon  the  dependent  variable  value. 
Dirichlet  boundary  conditions,  i.e.,  function  specified,  are  "strong" 
in  this  terminology  and  the  higher  the  order  of  the  derivative  boundary 
condition  applied  the  "weaker"  the  constraint  upon  the  dependent  varia- 
ble value  at  the  boundary.  The  "weak"  exit  boundary  conditions  referred 
to  here  are  in  effect  second  derivatives  set  to  zero  with  one-sided  dif- 
ferences. In  view  of  the  complexity  of  the  equation  system  involved  it 
is  not  a priori  clear  that  this  third  technique  of  linear  extrapolation 
would  workl  Certainly  all  the  suggestions  are  reasonable  and  the  choice 
is  simply  a question  of  computing  economies  weighed  against  level  of 
approximation  and  desired  accuracy. 
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This  third  technique  is  probably  the  least  inconvenient  to  apply  and 
will  be  used  unless  the  solutions  indicate  that  this  is  either  unstable 
or  inaccurate.  Availability  of  an  inviscid  code  to  compute  the  forebody 
flow  field  and/or  the  computational  costs  would  indicate  which  of  the  two 
other  strategies  would  be  adopted  should  the  linear  extrapolation  scheme 
be  unacceptable. 

The  upstream  boundary  conditions  are  available  from  conventional 
propulsion  and  interior  ballistics  codes;  see  for  example  Anderson  et  al  , 
->9  30  31 

Dahm  and  Anderson“‘ , Celmins  and  Gough'  . The  potential  inaccuracies 
of  these  one-dimensional  codes  is  a source  of  concern.  For  instance 
leakage  flow,  the  subject  of  the  detailed  analysis  under  consideration, 
must  obviously  have  an  effect  on  the  breech  pressure  level  and  hence  the 
projectile  velocity,  as  indeed  so  must  the  gap  friction.  Thus  there 
should  be  a coupling  between  the  local  gap  flow  analysis  and  the  overall 
interior  ballistics  predictions.  However,  at  least  in  the  initial  phases 
of  the  detailed  gap  analysis,  the  coupling  between  the  one-dimensional 
breech  analysis  and  the  gap  will  be  ignored. 

The  boundary  conditions  for  the  turbulence  model  will  be  treated  as 
follows.  At  the  wall  two  options  are  open.  In  the  first  the  grid  point 
definition  normal  to  the  wall  is  either  adequate  or  made  to  be  adequate 
to  define  the  viscous  sublayer  and  hence  boundary  conditions  of  zero 
turbulence  at  the  wall  are  physically  reasonable.  The  difficulty  with 
this  approach  is  that  apart  from  the  required  grid  point  definition,  the 
physics  of  low  Reynolds  number  turbulence  must  be  modeled  in  a reasonable 
manner  by  governing  turbulence  equations.  The  alternative  approach  is 
to  force  the  turbulence  equations  to  give  predicted  levels  at  the  first 
grid  point  away  from  the  wall  which  are  consistent,  in  some  sense,  with 
the  well-known  boundary  layer  of  the  wall.  The  law  of  the  wall  can  be 
written  for  the  axial  velocity  u for  instance  as 

u+  = f(y+) 

where 

+ \/2 
u = u/u  , y = yu  /v  u = (r  /p) 

T T L W 

and  y+  is  the  distance  normal  to  the  wall  and  xw  is  the  wall  shear  stress. 
The  function  f is  a known  universal  relationship  valid  for  smooth,  rough 
and  transpired  turbulent  boundary  layers.  Equivalent  relationships  are 
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available  for  temperature  species  and  turb  .ience  quantities.  The  law  of 
the  wall  can  be  applied  much  as  given  and  the  solution  of  the  finite  dif- 
ference equations  forced  to  satisfy  the  given  equation.  The  given  func- 
tional form  can  be  differentiated  and  derivative  conditions,  i.e., 
"weaker"  boundary  conditions,  used.  No  matter  how  implemented,  the  dif- 
ficulty with  this  whole  approach  is  that  the  validity  of  the  law  of  the 
wall  under  the  local  conditions  of  the  gun  barrel  gap  leakage  problem  is 
questionable.  The  law  of  the  wall  matching  technique  has  been  largely 
favored  in  the  past  because  it  economizes  on  grid  point  numbers  and  also 
because  presently  little  is  known  about  low  Reynolds  number  turbulence. 

Law  of  the  wall -based  constraints  for  turbulence  model  equations  have 
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been  developed  by  Buggeln  and  McDonald  . 'lhis  type  of  law  of  the  wall 
constraint  can  be  arranged  to  allow  the  effect  of  surface  roughness  to 
enter  and  it  will  be  used  in  the  numerical  model. 

The  downstream  boundary  conditions  for  the  turbulence  model  equations 
would,  if  possible,  be  of  the  very  weak  interior  extrapolation  type  re- 
ferred to  earlier.  As  before,  setting  the  second  derivative  of  the  de- 
pendent variable  to  zero  should  enforce  rather  weak  constraints  upon  the 
solution.  Should  this  turn  out  not  to  be  the  case  then  the  indications 
are  that  the  downstream  boundary  should  be  moved  closer  to  the  quiescent 
air  outside  the  weapon.  As  far  as  upstream  boundary  conditions  on  the 
turbulence  are  concerned,  little  information  is  available.  Here  all  that 
can  be  done  until  experimental  guidance  is  obtained  is  to  make  predictions 
for  various  plausible  upstream  turbulence  levels.  For  instance,  turbulent 
energies  varying  between  0 and  say,  100%,  of  the  bulk  gas  velocity  would 
seem  reasonable.  Further,  a turbulent  length  scale  varying  between  zero 
and  the  barrel  radius  seems  to  bound  the  possible  scales  involved.  With- 
in these  bounds  all  that  can  be  done  is  to  examine  the  sensitivity  of 
the  solutions,  in  particular  the  wall  heat  transfer  to  these  parametric 
variations.  In  view  of  the  high  acceleration  levels  near  the  gap  inlet 
it  is  suspected  that  the  gap  fluid  dynamics  will  be  relatively  insensi- 
tive to  the  breech  turbulence  levels. 

F.  Computer  Storage  and  Related  Problems 

In  developing  the  arguments  earlier  for  a block  implicit  approach 
it  has  been  claimed  and  assumed  that  the  resulting  linear  system  can  be 
solved  efficiently,  i.e.,  with  at  least  the  same  order  effort  as  other 
competitive  schemes,  and  that  no  unreasonably  large  storage  demands  are 
made.  Here  it  is  supposed  that  the  investigation  will  result  in  a com- 
puter program  which  will  run  on  a CDC  7600  generation  machine.  Prior 
20  21 

work  ’ indicates  that  an  alternating  direction  decompositon  (ADI)  or 
matrix  factorization  scheme  as  it  is  sometimes  called,  clearly  can  result 
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in  a very  efficient  technique  for  solving  the  linear  equation  system. 

Here  we  distinguish  between  ADI  and  matrix  factorization  terminology  by 
reserving  the  term  ADI  to  schemes  where  the  intermediate  level  solutions 
represent  consistent  approximations  to  the  governing  partial  differential 
equations.  In  this  way  physical  boundary  conditions  can  be  applied  at 
the  intermediate  levels,  a critical  factor  contributing  to  the  success 
of  the  schemes.  Matrix  factorization  includes  both  ADI  and  the  other 
schemes  which  do  not  possess  this  consistency  property.  Consistency  is 
felt  to  be  extremely  important  and  inconsistent  schemes  are  not  recom- 
mended here.  What  is  not  widely  recognized  is  that  this  ADI  decomposi- 
tion has  a very  profound  beneficial  effect  upon  the  storage  requirements. 
Indeed,  using  this  decomposition  usually  only  the  number  of  lines  in  any 
one  space  dimension  equal  to  the  bandwidth  of  the  spatial  difference 
molecule, together  with  the  number  of  time  levels  involved  in  the  temporal 
difference  molecule,  need  be  in  fast  core  memory  at  any  point.  As  a prac- 
tical matter  then  the  core  storage  requirements  can  be  kept  to  a minimum, 
making  extensive  use  of  peripheral  storages.  The  line  matrix  block  elimi- 
nation usually  provides  sufficient  time  for  buffered  input -output  trans- 
fers to  be  completed  prior  to  the  code  requiring  the  stored  arravs. 
Block  data  transfers  rather  than  point  by  point  transfer  can  and  should 
be  used  to  speed  up  the  input-output  process.  Thus  it  is  not  expected 
that  storage  will  be  a major  problem  if  an  ADI  scheme  is  adopted. 

G.  Available  Codes 

A number  of  documented  codes  exist  which  must  be  considered  as 
potential  candidates  for  the  calculation  of  the  leakage  problem.  In 
selecting  and  judging  them  there  are  a number  of  requirements  which  must 
be  met:  the  code  should  be  able  to  handle  viscous  turbulent  flow  and 
finite  rate  chemistry,  thus  resolving  the  reaction  zone,  and  consequently 
exhibit  accuracy  over  a range  of  Mach  numbers.  Desired  added  features 
? include,  but  are  not  limited  to,  multiphase  and  multispecies  capability, 

and  ability  to  handle  multidimensionality  of  the  flow  geometry.  It  is 
desirable  that  the  code  not  exceed  the  large  core  storage  capacity  of  the 
CDC  7600  at  BRL,  and  running  times  per  case  should  be  less  than  one  hour 
since  the  code  is  intended  for  parametric  studies.  The  algorithm  must 
have  been  checked  for  accuracy  and  stability  under  typical  running  con- 
ditions. Finally  the  code  should  be  able  to  handle  efficiently  and  eco- 
nomically stiff  systems  arising  out  of  flows  with  realistic  chemistry. 


There  are  four  broad  categories  of  codes,  documented  and  generally 
available,  which  have  to  be  considered  for  the  reactive  leakage  flow 
problem:  the  Los  Alamos  family  of  codes  based  on  the  ICE  and  ALE 
techniques,  those  derived  from  variants  of  MacCormack's  scheme,  Spalding's 
effort  and  MINT  developed  by  Briley  and  McDonald. 
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Butler  has  reported  on  calculations  with  different  versions  of 
ICE.  Satisfactory  results  have  also  been  claimed  for  treating  a bank  of 
gas  dynamic  lasers,  flame  propagation  and  the  reactive  flow  inside  a com- 
bustion engine.  This  method  is  applicable  to  the  full  two-dimensional, 
time-dependent  Navier-Stokes  equation  with  species  transport,  mixing  and 
chemical  reaction  between  the  constituents.  The  RICE  code,  a reactive 
version  of  ICE,  used  in  these  calculatios,  is  based  on  a semi-implicit 
technique  which,  by  modifying  the  time  step  restriction,  permits  a more 
efficient  calculation  to  be  performed.  The  code  and  the  calculational 
results  can  be  criticized  on  a number  of  accounts,  and  in  the  present 
problem  the  principal  problem  would  most  likely  be  the  undue  stability 
restrictions  on  the  time  step.  As  a general  observation  with  scant  ex- 
perimental data  available  for  comparison,  it  is  rather  difficult  to  pass 
judgment  on  the  physical  accuracy  of  the  reported  results.  Flow  gradi- 
ents are  treated  by  numerical  approximation  and  if  this  ij  not  done  suf- 
ficiently accurately, either  as  a result  of  poor  choice  of  assumed  func- 
tional form  of  the  solution  or  by  an  inadequate  mesh,  then  poor  or  false 
predictions  of  the  physical  process  can  result.  It  is  felt  that  exten- 
sive validation  with  respect  to  measurements  need  to  be  performed  before 
any  of  these  codes  can  be  considered  for  alternate  applications. 

MacCormack's  scheme  has  been  shown  to  yield  acceptable  results  for 
reactive  flows  for  a number  of  different  geometries.  The  most  efficient 
version  of  this  explicit  method  uses  a formulation  whereby  the  chemistry 
is  split  off  from  the  hydrodynamic  part  of  the  problem.  Efficiency,  accu- 
racy and  economy  of  storage  have  been  claimed  for  this  technique  and  demon- 
strated on  a number  of  practical  problems.  MacCormack's  approach  runs 
into  a number  of  difficulties  and  restrictions  when  applied  to  flow  geom- 
etries where  the  nature  of  the  flow  undergoes  a drastic  change,  such  as 
a change  from  hyperbolic  to  parabolic  type  within  a boundary  layer. 
Running  times  for  the  problem  of  interest  could  be  quite  long  due  to  the 
time  step,  stability  restriction.  Also,  the  technique,  when  tried  on  a 
two-phase  problem, has  not  led  to  completely  satisfactory  results. 

Spalding  and  his  coworkers  have  been  responsible  for  the  creation 
of  a number  of  reactive  hydrocodes,  some  of  which  include  the  effect  of 
turbulence  on  the  reactive  flow.  However,  from  the  point  of  view  of  the 
leakage  problem,  these  codes  have  limited  applicability  since  they  pri- 
marily address  the  problem  of  steady  flow  at  low  Mach  numbers.  One  time- 
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suited  for  the  present  problem  is  not  generally  available.  Even  if  this 
code  were  available  the  pressure  treatment  is  somewhat  awkward  and  approxi- 
mate, making  it  an  unlikely  choice  for  the  present  problem. 

The  MINT  code  is  based  on  a fully  implicit  technique  and  is  a three- 
dimensional,  time-dependent,  compressible  multicomponent  code  which  has 
been  exercised  on  a variety  of  combustor  problems  with  satisfactory 
21 

results  . It  fulfills  the  criteria  set  forth  at  the  beginning  of  this 
section  and  it  can  be  adapted  for  the  purpose  of  this  analysis.  The 
modifications  needed  include  an  adaptive  coordinate  system,  a two  equa- 
tion turbulence  model  and  characterization  of  the  burning  propellant  bed. 

Perhaps  more  importantly,  present  familiarity  with  this  code  and  our 
confidence  in  its  reliability  of  accurately  modeling  some  aspects  of  two- 
phase  flow  phenomena  suggest  that  MINT  should  be  the  departure  point  in 
the  present  modeling  effort.  It  is  anticipated  that  compression  wave 
definition  may  pose  some  problems,  but  at  this  juncture  it  is  felt  that 
this  shortcoming  is  outweighed  by  the  advantages  that  this  approach  offers. 

Several  other  developmental  efforts  should  be  mentioned  at  this 

juncture.  Gough^5  is  in  the  process  of  extending  his  one-dimensional 
two  phase  flow  of  the  gun  tube  to  two  dimensions  using  the  Lax-Wendroff 
algorithm.  Experimental  verification  of  his,  as  well  as  all  other  codes, 

needs  to  be  performed.  Finally,  the  Calspan  model  , which  computes 
rather  than  takes  a specified  ignitor  input,  has  been  used  to  simulate 
a 155  mm  howitzer  with  a modicum  of  success.  However,  like  Gough's 
model  it  underpredicts  the  breech  pressure  as  well  as  the  muzzle  velocity, 
leading  one  to  question  the  physics  as  well  as  the  mathematics  of  the  code. 

H.  Experimental  Validation 

Little  or  no  existing  experimental  evidence  appears  available 
for  use  in  validating  the  code.  This  state  of  affairs  must  be  rectified 
before  confidence  in  the  analysis  can  be  obtained.  Nonreacting  experi- 
mental studies  can  be  performed  relatively  cheap  if  the  flow  is  steady. 
While  not  validating  the  transient  capability  of  the  analysis, such  tests 
do  prove  a very  necessary  check  since  there  would  be  little  hope  of 
treating  the  transient  problem  if  the  steady  flow  could  not  be  accurately 
predicted.  An  insulated  wall  with  a moderately  heated  gas  would  enable 
wall  temperature  levels  to  be  predicted  and  compared  to  data  for  various 
gap  heights  and  Reynolds  numbers.  With  high  response  instrumentation 
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the  wall  or  projectile  could  be  translated  to  introduce  some  transients 
into  the  flow.  Once  again  wall  temperature  would  be  a stringent  test  of 
predictive  accuracy. 


Reactive  flow  tests  would  be  extremely  difficult  and  probably  wall 
temperature  measurements  at  fixed  locations  on  the  bore  of  a weapon  would 
be  the  most  one  could  expect.  Even  this  very  modest  type  data  would  be 
difficult  to  obtain  but  very  valuable  in  view  of  the  realistic  environment. 

In  both  the  hot  and  cold,  steady  and  unsteady  and  experimental 
studies,  mapping  of  the  flow  field  would  be  highly  desirable.  Realistic- 
ally speaking,  however,  such  measurements  are  probably  not  feasible  with- 
in the  constraints  of  time  and  effort  required  to  obtain  them  at  present. 

V.  Recommendations  for  Future  Work 

Our  principal  recommendation  is  to  adopt  the  MINT  code  to  the  needs 
of  the  interior  ballistics  modeling  problem  subject  to  the  extensions 
and  restrictions  discussed  in  this  report. 

Future  generalizations  follow  from  our  discussion.  These  include, 
but  are  not  limited  to  an  accurate  description  of  the  effect  of  the  par- 
ticles on  the  flow  as  well  as  consideration  of  particle-particle 
interactions.  Realistic  chemistry  to  take  cognizance  of  the  propellant 
bed  dynamics  as  well  as  the  interplay  between  burning  rate  and  turbulence 
will  have  to  become  an  integral  part  of  the  model. 

The  boundary  treatment  should  reflect  the  heat  losses  through  the 
wall  which  may  influence  the  flow  definition.  With  the  availability  of 
sufficient  computer  storage,  three-dimensional  effects  such  as  balloting 
and  swirling  flow  induced  by  the  rifling  of  the  gun  tube  should  and  can 
become  an  integral  part  of  the  model. 


Figure  1.  Schematic  of  Projectile  in  a Smooth  Bore  Gun  Barrel  at 
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